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TECHNICAL PUBLICATION 


A REVIEW OF THE GINZBURG-SYROVATSKII’S GALACTIC COSMIC-RAY 
PROPAGATION MODEL AND ITS LEAKY-BOX LIMIT 

1. INTRODUCTION 


Galactic cosmic rays (GCRs), while mostly protons, are energetic nuclei representing per- 
haps all of the nuclides of the periodic table. Their origin remains unclear even as the centennial of 
their discovery by Victor Hess is celebrated in 2012. In his pioneering balloon-borne experiments, 
Hess discovered that the level of ionizing radiation was actually higher up in the atmosphere than 
at sea level (e.g., by a factor of 2 at a height of 5 km). He concluded that the increase was due 
to ‘radiation’ penetrating the atmosphere. Robert Millikan confirmed Hess’s discovery in 1925 and 
coined the term ‘cosmic rays’ (CRs). Hess was awarded the Noble prize in physics in 1936 for his 
discovery of ‘cosmic radiation.’ (Hess shared the 1936 prize with Carl Anderson for Anderson’s 
discovery of the positron.) The histories and disciplines of nuclear and particle physics have been 
interwined with CR physics ever since. The 1947 discovery of the pi-meson in ‘atmospheric CRs’ is 
a telling example of this common heritage. 

Today, a host of balloon-borne, space-borne, as well as ground-based observations of (pri- 
mary) CRs reveal they are nuclei of regular (i.e., baryonic) matter from the lightest, hydrogen, through 
the very heavy elements like uranium. In addition, CRs also appear to include electrons, positrons, 
and antiprotons. Their arrival, for the most part, i.e., better than 1 in 10 4 , is isotropic. The chemical 
abundances — elemental as well as isotopic — of these primary CRs appear to reflect the chemical 
composition of massive stars’ winds (as opposed to solar system-like abundances, for example). 

Their arrival spectrum near Earth’s orbit reveals an almost structureless power law spanning 
many decades of energy, down to ~10 18 eV. Primary CRs arriving at the boundaries of the helio- 
sphere with energy below a few hundred MeV/nucleon are unable to penetrate the outward-flowing 
solar wind plasma. Those that do lose a fraction of their energy as they drift and diffuse into the 
solar wind — an effect known as solar modulation. Near the inner heliosphere and Earth’s orbit, the 
primary CR nuclei arriving at the top of the atmosphere interact with its atomic constituents (mostly 
carbon, nitrogen, and oxygen), losing (as they cascade through), in most cases, all of their kinetic 
energy while producing a host of secondary radiation and particles (like gamma rays, neutrons, and 
pions). Very few primary CRs make it to sea level. Their secondary products, however, like muons, 
are able to penetrate the hard surface of the Earth down to depths of tens to hundreds of meters. 


From an astrophysical perspective, some of the outstanding fundamental questions regard- 
ing CRs pertain to the nature and location of their source(s) as well as to the precise mechanism(s) 
of their acceleration. Their journey after synthesis and acceleration — their propagation in the inter- 
stellar medium (ISM) — appears to be reasonably well understood. Currently, the so-called diffu- 
sive shock acceleration theory appears to be the accepted physical theory behind their acceleration. 
The ‘standard’ source is usually taken to be a supernova explosion. However, galactic superbubbles, 
where a collection of massive OB stars tend to congregate in close proximity in space and time, are 
increasingly becoming a more likely source. 

The propagation stage of that cosmic journey will be the focus in this Technical Publication. 
The medium of propagation is the ISM and, for purposes here, this is idealized as a gas of mostly 
hydrogen and helium (He) with a number density on the order of 1 per cm 3 and a kinetic tempera- 
ture of ~10 4 K, permeated by a weak (~pG) but turbulent magnetic held. The CRs themselves are 
expected to contribute to the energy density of the ISM (which is ~2 eV/cm 3 ), but these and similar 
effects will be ignored here. Finally, the particles’ transversal of this medium is describe only statisti- 
cally, and classically, taking into account only basic energy gain and loss processes as well as regular 
(baryonic) particles’ sinks and sources. 
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2. THE GINZBURG-SYROYATSKII EQUATION 


In transport and acceleration of energetic (suprathermal) charged particles in space and 
astrophysical plasmas, the relevant length scales are typically much larger than the particle’s Larmor 
radius. As a result, the motion of these particles is effectively determined by the structure of the elec- 
tromagnetic held via both its regular and stochastic components. In the absence of particle-particle 
collisions (the low density limit), the electromagnetic held is also responsible for accelerating the par- 
ticles up to relativistic energies. The huctuating component of the held can scatter particles, especially 
when the particle’s Larmor radius is comparable to the wavelength of the scattering hydromagnetic 
wave, ‘resonant scattering.’ If scattering is frequent and strong, it can isotropize the particle’s density 
function and, in the diffusive limit, the motion of those particles can then be described statistically. 

As CR nuclei traverse the ISM and interact with its constituents and helds, they can lose or 
gain energy as well as suffer nuclear interactions. The Ginzburg-Syrovatskii (G-S) equation 1 follows 
the evolution of a CR particle density distribution function, N(E, r,t ), that evolves in space and time 
once created at t - 0 and by some source somewhere in the galaxy. N is assumed to depend explicitly 
on the kinetic energy, E, (or momentum) of the CR particle. For a particular CR nuclide i, the G-S 
equation for N t (E, r,t ) is a continuity equation in phase space and reads (eq. (14.8) in ref. 1): 


dN i 

~dt~ 


v-(A VJV)+^(W) 


1 


2 dE' 


{Wi) = Qi 


-Pi N i + Y,Pi N k- 

k<i 


( 1 ) 


On the left-hand side of the G-S equation, ZX is the diffusion coefficient (or tensor) for species 
i in the ISM, which, in general, can depend on E, t, and f. The coefficient or function b i = dEUlt 
descri bes the systematic energy change, A E, while d t describes the fluctuations in this change, i.e., d t (E) 
=d/dt(AE) 2 . On the right-hand side of the G-S equation, Q t (E,r,t ) is the source intensity, / , of the 
CR nuclide i, where, for an isotropic distribution of particles, N (>E) = 4 n\ IAE) / v dE, where v is 
the particle’s speed. p t is the probability (per unit time) for a CR nuclide i to suffer a collision in the 
ISM with a cross section o t , i.e., p(= nva h where n is the ISM number density. In the last term on the 
right-hand side, p,' is the probability per unit time for the appearance of a nuclide i due to collisions 
of nuclide k in the ISM, and hence, p- can be written as = nvo ik , where a ik , is, for example, the pro- 
duction cross section for nuclide k from spallation reactions of nuclide i. In this notation, index i- 1 
refers to the nuclide with the highest atomic number. Note that the G-S equation also applies to CR 
electrons with some modifications. 1 


In eq. (1), terms related to further addition and subtraction of species can be added as well 
as terms related to further energy losses and gains. But note the clear distinction of space-time 
evolution, ( dN/dt ) and (V *(DVN), from chemical transformation in which energy appears only as 
a parameter rather than as an independent variable, if the energy loss or gain terms are ignored. 
Actually, for CR particles with kinetic energy greater than a few GeV/nucleon, energy loss (or gain) 
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mechanisms can, for propagation purposes, be neglected to a good approximation. As will be seen 
shortly, the transformation part can be separated and rewritten in terms of the amount of ISM mat- 
ter, x, traversed by the CR nucleus as 


dNf 

dx 


-Xi 


k<i 


( 2 ) 


where, under the homogeneity assumption, v = nvt which is typically expressed in units of g/cm 2 , and 
n being the (average, in this case) matter density of the ISM. In equation (2), X t is the nuclear inter- 
action mean free path for species i in the ISM and X lk in the interaction mean free path to produce 
species k upon the interactions of species i with the ISM; X i and X ik are also typically expressed in 
g/cm 2 . Note that X°c 1/cr. Also implicity assumed in equation (2) is that N-(x = 0) = q i , qj being the 
fraction of nuclide i in the CR source and that N- is only a function of x. 

Note that eq. (2) is a special case of eq. (1) for Z) =b i =d i = 0 and for Qj-Qj(t)= qfiit). Essen- 
tially then, eq. (1) follows the evolution of N t as it evolves in time (or x) and diffuses in space while 
suffering, while not necessarily in parallel; however, nuclear transformations due to interactions in 
the ISM which can alter the particular species N t describes (for a thorough treatment of nuclear 
spallation processes of CRs, see ref. 2.) The above space-time versus chemical transformation separa- 
tion (which will allow for a general solution to the G-S equation to be written down) is made possible 
by a number of assumptions, most important of which is the assumption that the space-time and 
transport parameters D i , b k d { do not depend on nuclear species while the transformation param- 
eters Pj and p- do not depend on either time or space coordinates. Before the general solution and 
various simplifying assumptions of the G-S equation are examined, it is worthwhile to appreciate, at 
least qualitatively, the justification as to what makes a diffusion-like equation appropriate to model 
the propagation of CRs in the ISM. 
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3. THE DIFFUSION APPROXIMATION 


Phenomenologically, the propagation of GCRs in the ISM is thought of as an electrodynamic 
interaction of a relativistic charged particle in a turbulent, magnetized cosmic plasma. The subject, 
of course, is rich in history and accomplishments, as well as in complexity. It is typically formu- 
lated within the context of a wave-particle interaction using the tools of weak turbulence theory. On 
a lower level still, i.e., more formal, one can think of it as a problem (and a solution) belonging to 
time-dependent perturbation theory. None of these tools or theorems will be recited here. Rather, the 
aim is to gain insight into how one can abridge a notion about the motion of a fast charged particle 
in a turbulent field with the diffusion approximation so as to arrive at the G-S differential equation, 
i.e., making the passage from the basic notion of motion in a turbulent magnetic field to a governing 
equation (like the G-S equation) with reductionist and predictive powers. 

Diffusive motion due to some randomness in the problem, e.g., Brownian motion, is a well 
studied notion. 3 At some level, one envisages a particle making a small and discrete number of ran- 
dom walks, due to irregular scattering, for example. A phenomenological description amounts to 
writing down a probability density function that can give a measure of where the particle would be 
and with what probability. A Bernoulli distribution function, for example, for the particle’s motion 
in one-dimension and for a small number of discrete steps can deliver that. (Keep in mind that the 
Bernoulli probability density function (PDF) has an asymptotic Gaussian form in the limit of a 
large number of steps. This is helpful because the Gaussian form also happens to be a solution to 
the simple diffusion equation.) When one proceeds to formulate the problem for a large number of 
random walks, the process (solution) begins to take on a simple form. Further, instead of addressing 
the motion of a single particle and its PDF, one thinks of a distribution function of a large number 
of particles that share the same initial conditions. To make the passage then to a differential equation 
governing the motion (with the appropriate initial and/or boundary conditions), one assumes some 
time step At that is large enough so that the particle will have time to make a large number of steps 
but short enough so that the rms of the displacement is small compared with large-scale motion of 
the particle. Now, in the limit At^O and by taking a Taylor-series expansion of our macroscopic 
distribution function N around At, one can arrive 3 at a differential equation for N that resembles 
a^diffusion equation, i.e., a partial differential equation (PDE) first order in time and second order in 
space. 


One can also arrive at a similar PDE by invoking the macroscopic Fick’s law along with the 
continuity equation. Fick’s law reads: 


j = -mN, 


(3) 


where J is flux and D is the diffusion coefficient (or tensor). The continuity equation for N reads 


dN 

dt 


+ V/ = 0. 


(4) 
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Now, combining eq. (3) with eq. (4) for N yields: 


^-h-V-(Z)VAT) = 0. (5) 

ot 

One can add a convection term, a source term, energy loss/gain terms, particle loss/gain terms, etc. to 
eq. (5) to arrive at a PDE more or less resembling the G-S equation. Note that eq. (5) is a linear PDE 
for N. (The same linear PDE (for N above), however, can be nonlinear if formulated in the ‘wrong’ 
variable, e.g., for v where J = v N, 4 as in eq. (28).) 

A Langevin 5 ’ 6 approach, in contrast, would specify instead, a stochastic differential equation 
(SDE) for the particle’s own random position, x, for example, as: 5 ’ 6 


dx-C(x,t) dt + jD(x,t) dW , (6) 

where for the scalar random variable x{t), C{x,t ) is its drift (or mean) term and D(x,t) is its diffu- 
sive (or dispersive) term, and dW is a Wiener process (i.e., a Gaussian white noise). While the two 
approaches are numerically different, they are nonetheless considered equivalent. 5 
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4. PARTICLE- WAVE INTERACTIONS 


Being very qualitative here, a few words about the so-called quasi-linear approximation is 
in order. When a group of particles interact with a wave (magnetic held) and they remain in phase 
with the wave over many cycles, the interaction is termed a resonant wave-particle process. In this 
process, the particles do have time to exchange energy with the wave while the distribution function 
will only suffer slow changes. Such changes are attributed to a diffusion process termed quasi-linear 
diffusion. The process is kinetic in nature, with the characteristic plateau and tail formations in the 
distribution function. Quasi-linear theory (QLT) has been most successfully — but still not without 
some difficulties — applied to the transport of solar energetic particles (SEPs) as well as CR particles 
in the interplanetary magnetic held (IMF ). 7 ’ 8 

In the case of the IMF, the magnetic held is modeled as having a uniform component B () (an 
Archimedian spiral emanating from the Sun) and a fluctuating irregular component SB. The motion 
of a typical SEP or CR particle then in this held will comprise of a smooth, regular component along 
B 0 and a chaotic component due to SB. The latter component results in the scattering of the particle 
in pitch-angle p (p is dehned as the cosine of the angle between the particle velocity vector and the 
magnetic line of force). If the scattering is assumed to be ‘sufficiently strong,’ then a diffusion process 
is a good approximation for the evolution of N. To quantitatively describe this diffusion process, 
a Fokker-Planck equation can be constructed that follows the evolution of N. 

Determination of the diffusion coefficient from knowledge of the magnetic huctuations is the 
essence of any statistical description of the transport process. In QFT formalism! s), the assumption 
is that those huctuations affect the motion of the particle on timescales slower than the characteristic 
timescale of the motion in the ambient held in general. In other words, a perturbative approach is 
assumed to hold true and the motion can be looked at as being uniform to the zeroth approxima- 
tion plus an irregular one to the hrst approximation. In the zeroth approximation, the particle ‘sees’ 
a stationary magnetic held, so a magnetostatic approximation can be used to derive the diffusion 
coefficient in terms of the observable statistical properties of the IMF like its mean value and power 
spectrum. This is an example of a Kolmogorov formalism wherein an incomplete knowledge of the 
nature and geometry of the magnetic huctuations is replaced by a hierarchy of correlation functions 
and their Fourier transforms. 

Other processes that can occur, but typically not included in (early) standard QFT of CR 
transport, are nonlinear wave- wave interactions and wave-particle-wave interactions . 9 The former, 
also a resonant process but no particle resonances are involved, can give rise to the so-called para- 
metric instabilities, where two waves beat together to interact with a third wave with the interaction 
becoming unstable above a certain threshold. The latter is both a resonant and a kinetic process 
where the particles keep in constant phase relative to the beat frequency of two waves, giving rise to 
the so-called plasma wave echoes . 10 
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5. GREEN’S FUNCTIONS SOLUTION TO THE GINZBURG-SYROVATSKII EQUATION 


The Green’s functions approach to a general solution of the G-S equation is appropriate here 
by virtue of the fact that the G-S equation, eq. (1), can be shown to satisfy the nonhomogeneous 
Sturm-Liouville equation which, in one independent variable x, reads: 

0[z{x)\ + f{x) = 0, (7) 

where O is a differential operator of the form 


0 [-] = 


d_ 

dx 



[■••]+ 2(X) , 


( 8 ) 


and pipe) and q(x) are some functions of x. Differential operators satisfying eq. (8) are called self- 
adjoint operators. One then seeks a solution, z(x), satisfying eq. (7) in addition to certain boundary 
conditions at endpoints a and b over the interval [a, b], with the help of another function (Green’s 
function). The sought-after Green’s function, G(x), is required to satisfy the homogeneous Sturm- 
Liouville equation: 

O[0(x)] = 0 (9) 

in addition to other conditions. 11 ’ 12 Most notable is that Q has to satisfy the same boundary condi- 
tions that z(x) satisfies and that G be a continuous function of x while its first derivative with respect 
to x be discontinuous. In terms of Q, the general (Green’s functions) solution to eq. (7) is now written 

z(x)= [ dx'G(x,x')f(x'), (10) 

where the integrated-over variable x is not exactly a dummy variable; rather, it is the variable that 
G' needs to be discontinuous in when p(x), eq. (8), is evaluated at. 

Back to the G-S equation, a general, Green’s function solution 13 for N satisfying eq. (1) can 
now be written as 


N i F ) = J 0 dx F ’ X ^ N i ( A ‘ ) ’ C 1 !) 

where x is an independent variable. Note, however, that x and t will no longer be independent under 
the homogeneity assumption and the N- (x) are functions of x only that satisfy 


dN- 

dx 


-°iN‘i + X °i j N S j + c i 8(x) . 
j<i 


( 12 ) 
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The source term is now a Dirac delta function in x with strength c t (to simulate relative abundance of 
the different species /' at synthesis). Collecting the interaction terms together and ignoring the energy- 
loss term, Q is now required to satisfy 


d_ 

dt 


V- (Z)V) + nv 


dx 


G(t,r,x) = c5(x ) , 


(13) 


in addition to the boundary conditions that N satisfies. Solving for N amounts to constructing Q, 
then solving eq. (11). Intuitively, one can think of Q as that key part of the diffusion model that 
embodies all space-time as well as x information since eq. (12) carries no space-time information, i.e., 
the standard notion of a Green’s function as an ‘influence’ function is clear here as well. 


Constructing a general (or a class of) Q, however, is not an easy task, and the carrying of 
space-time information will, to a large degree, be model dependent. As a prime example, in the lit- 
erature, 13 ’ 14 a steady-state model, i.e., dN/dt = 0, with cylindrical symmetry has been put forward, 
meaning a general Q is written down with space information as well as other astrophysical input, e.g., 
relative thickness and density of the galactic disc to halo, etc. The model assumes that CR particles 
can freely escape the galactic volume (disc + halo) by crossing the surface S that encloses this volume, 
i.e., the assumed boundary condition at the surface is that | 5 = 0 (as opposed to a reflecting bound- 
ary, dN/dz \ s = 0, where z is some space coordinate). Also, it is further assumed that Nj is continuous 
across the boundary of the disc and halo. While the complete formulation of this construction of Q 
is rather involved, for purposes here, the limit of the formulation in one-dimension is what is more 
instructive. It is helpful to think of that limit as when the diffusion process takes place mainly along 
a tube of magnetic field lines, and the only relevant space coordinate in this case will be the one along 
the field lines of the galactic magnetic field. Also, for purposes here, z (say z = 0) can be pinned down 
and the essential information will still be intact. In this case, and for a constant diffusion coefficient, 
D. Q takes the simple form: 


s,.(-- = o)~d 


1 + nvAO: ID 

& 1 


(14) 


Above, A is an 
cross section for 
<*i « D / ( n^vhj, ) 


‘area’ of the order of h g h h , disc (gas) and 
nuclide i, and n g is the disc’s matter density. 


halo heights, <r ; is the total reaction 
The above approximation is valid for 
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6. TRYING OUT DIFFERENT <?S 


One can construct a global model for which one derives the appropriate Q, as in the 
above model of Ginzburg and Ptuskin, 14 or one can try out various forms for Q and pick the 
one that seems to fit the modeling purpose. It is instructive to go through some examples, 
however, as they will help us connect with the various notions, models, etc. that are widely used — 
especially in data analyses. 


6.1 The Slab Model 

The simplest form one can pick is that of a Dirac delta function: 

£->(5(x-x 0 ). (15) 

This particular choice of Q carries no explicit space or time information. The only information it car- 
ries is that of a path length x = x 0 , i.e., species i is propagated through ISM matter equivalent to path 
length Xq. For this Q, eq. (1 1) becomes 


Ni(t,r)^N-(x 0 ). (16) 

N- (a' 0 j , in turn, is the solution to eq. (12) evaluated at x = x 0 . Incidentally, eq. (12) has a general 
solution of the form: 15 


N i (•*) = X[ ex p( +M 0] - N j (* = °) » 


(17) 


where M x is a species transformation matrix, i.e., nuclear chemistry information. Thus, with the 
choice (limit) of a Dirac delta function for Q, the solution is reduced to the so-called ‘slab’ model of 
GCR propagation. Also obvious is that this limit suppresses all spatial information about the diffu- 
sion process, i.e., nothing is being diffused, only species transformed. Or, diffusion was so rapid, sug- 
gesting space homogeneity throughout the transit process. The limit, however, does serve the purpose 
of pointing out the complete separability of terms carrying the space-time information (diffusion), 
Q, and those that are responsible for species transformation, N- (x) , in eq. (1 1). This feature, which 
comes about from the form of the G-S equation, is useful for studies of the GCR source composi- 
tion as it allows one to ‘propagate back’ to the source, insofar as species transformation is concerned, 

simply by solving for N- (x = 0) given N- (x () j : 


A', s (x = 0) = £[exp(-M,) 
j 



(18) 
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But since the diffusion process is irreversible (initial conditions about space-time are lost in the pro- 
cess), propagating back cannot refer to this part of eq. (11). Fortunately, the separability feature 
alluded to above and the fact that the various nuclear information (cross sections) are not expected 
to be space-time dependent, one can, with some confidence, speak of propagating back to the GCR 
source. 


6.2 The Venerable Leaky-Box Model 

Starting here with the results, then the steps will be justified. The Q in the leaky-box limit is 
that of a simple exponential: 

Q(x) = exp(-x / v esc ) , (19) 

where x esc is the so-called escape path length related to the escape time (from the galaxy) of GCR by 
a; = nvT , where n is average density. Back to eq. (1 1), with the above Q, one has a solution for 
N in the so-called leaky-box limit. One recognizes that, instead of a simple slab model, i.e., a singular 
v, here is a distribution of x (path length distribution or PLD). This is termed the ‘weighted slab’ 
model due to the appearance of Q in eq. (1 1) as a weighing function in x. Thus, the leaky-box limit 
(or model now) and the weighted slab model are mathematically equivalent in the absence of any 
energy losses, as assumed from the outset. 15 ’ 16 In the presence of such losses, one can still treat them 
equally (at least for numerical purposes) as long as energy appears as a parameter in the equation to 
be solved. 

Let us begin our justification by making the leaky-box ansatz that GCRs, for the most part, 
are confined to the galactic volume except for those GCRs that are endowed with enough kinetic 
energy to escape the confining effect of the galactic magnetic field (see fig. 1). Note also that the char- 
acteristic size of our galaxy is of the order of 100 kly (kilo light years). Since these particles travel 
at almost the speed of light, the characteristic transit time should then be of the order of thousands 
of years. Cosmic-ray data, however, as inferred mostly from the secondary to primary ratios which 
are directly affected by the amount of ISM matter traversed (of the order of 10 g/cm 2 ) suggest that 
this timescale is actually of the order of 10 My. Hence, CRs travel in long-winded (diffused) paths 
between their source(s) and detection. Also note that CRs with energies above 10 19 eV have a gyro- 
radius in the ISM’s magnetic field (strength on the order of pG) larger than the galactic diameter 
of -100 kly; they can escape. Those that do escape are assumed to be able to do so at any point on 
surface S enclosing the galactic volume. 

There is a subtle point here with regard to boundary conditions. In the diffusion model, the 
boundary condition is such that N drops to zero at the boundary point (Dirichlet boundary condi- 
tion), any point on S. Now, in this leaky-box ansatz, a strong reflection, i.e., the gradient of N drops 
to zero at any point on S (Neumann boundary condition), is assumed for most of the CR particles 
except for those that can escape, of course. As such, i.e., different boundary conditions, the diffu- 
sion model and the leaky-box model cannot, strictly speaking, be considered related (even though a 
fundamental relationship in what follows will be established). What this means is that the leaky-box 
limit can be considered a truly mathematical limit (however physical or unphysical that limit may be) 
of the diffusion model at any point inside S; but not on S itself. (It is not clear what implications, if 
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Figure 1 . A schematic (not to scale) showing the long-winded path of a high-energy 
cosmic-ray particle in the Milky Way as it makes its way out of the ‘box.’ 


any, this ‘technical’ point may have on conclusions about the size of the confining volume reached 
via the leaky-box model.) Recall that in addition to this, in the leaky-box model, one is also assuming 
spatial homogeneity. 

Now, the first link that ought to be established between the two models is the simple exponen- 
tial form of Q (in the diffusion model), or equivalently, the simple exponential PLD characteristic 
of the (standard) leaky-box model. One arrives at the leaky-box equation by replacing the diffusion 
term with an escape term (the homogeneity assumption), 


-v (i>Viv)^+^ , 

^esc 

and with no explicit time dependence, dN / dt - 0 , eq. (13) for N will read, 


( 20 ) 



dN 


+ nv — — 
ox 


q , 


( 21 ) 


where both the density and source terms now denote average density and average source so as to 
be consistent with the homogeneity assumption. Equation (21) is the leaky-box equation without 
energy loss/gain or decay terms. Now, Q must satisfy a similar equation: 


T 


esc 


+ nv — 
dx 


■ S(x) . 


( 22 ) 


where, for simplicity, the strength of the delta function source term has been set to unity. Solving 
eq. (22) for Q gives, 
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(23) 


Q{x) = exp(-x / a 6SC ) , 

which, after inserting it back into eq. (11), yields the solution of the G-S equation in the leaky-box 
limit. The escape path length x esc will be related to the diffusion coefficient as will be shown when the 
connection between the two models is established more firmly. 

But, before that is done, the two underlying assumptions — steady state and homogeneity — 
need some attention. Steady state is readily argued for on the basis of relevant GCR propagation times- 
cales; i.e., relevant times are much longer than it takes a typical GCR particle to traverse the confining 
volume once, implying that a steady-state situation can be safely assumed. Homogeneity deserves 
a little algebra, and will be looked at from two different angles. 

Species transformation aside, recall the separability feature of the G-S equation and examine 
the following term: 

-V- (.DViV(F)) = q , (24) 


where now N is a function of space only due to the steady-state assumption. Now, for a homoge- 
neous volume V' (enclosed by S' rather than S because of the different boundary conditions), and 
after averaging over that volume, gives for the left-hand side of eq. (24): 


-V- (f)VN(r )) = J </V[-V- (DVN)] / V' 
= jdS'(-DVN)n/V' 

= J<«'A scK. 


(25) 


where n is a unit vector _L S ' . Integrating over S' the J’s vanish except for ./ esc , which is the leaky- 
box ansatz. Now, if J esc = v esc N and also assuming that where the particle escapes (on S') is inde- 
pendent of S ' , then 


-V-(i)VATr)) 




S K 


esc 



(26) 


In eq. (26), the volume V’ ~ S' r' sc , where r e ' sc is the characteristic escape distance and T esc is 
the characteristic escape time. Averaging over V' on the right-hand side of eq. (24) gives q . Thus, 
eq. (21), after inserting back the species transformation term, is derived under the homogeneity 
assumption. 
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Also, for a constant diffusion coefficient, D, the diffusion term in the G-S equation under the 
leaky-box ansatz reduces to A/const where the constant (with the dimension of time) can be identi- 
fied with T esc . Start by writing 


-V- (DVN) | esc = V- Av esc , (27) 

which is valid for any point either on or inside the enclosing surface, be it S or S' ■ The right-hand 
side of eq. (27) can be expanded to 


V ■ Nv = NV ■ v + v ■ V N 

esc esc esc 




v — V 
esc esc t 


N . 
= — O 

D 


(28) 


where |V- v esc - v 2 sc j = O, which has the look of a ‘potential.’ Rewriting the diffusion equation in 

terms of v and at v gives 
esc 


dv 


esc 


dt 


= v(V-v -v 2 ) 

\ esc esc / 


= VO. 


(29) 


Now, if escape is not an explicit function of time, i.e., 3v esc / dt = 0, then O = const, which for a con- 
stant D gives, 


-V • (DVN) | esc = N / const , 


(30) 


which is again the leaky-box ansatz. 

The difference between the two derivations, using S' and O (even though they both lead to 
the same leaky-box limit) is that the question of boundary conditions does not come up when using 
O. What this means, though, is that the identification of the constant in eq. (30) as that related to 
some (constant) characteristic escape time seems valid, while not unique. This points again to the 
question of the different boundary conditions in the two models and why one should be somewhat 
careful when using conclusions of one model to infer or support assumptions about another. 

6.3 Nonstandard Leaky-Box Q 

As seen above, the Green’s function for the standard leaky-box limit is that of a simple expo- 
nential, eq. (19). Variants 17-19 of this simple function have also been used and found to fit some of 
the observed CR data. These include the double exponential, Q- where both Q x and Q 2 have 
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the form of eq. (19), but with different x esc . This is the so-called nested leaky-box model and the 
physics behind it is that of a ‘shrouded’ CR source. Also, there is the truncated exponential PLD: 


G(x) = < 


0 , 

exp[-(x-<5/x esc ] 


x < 8 
x>8 


(31) 


where one assumes a significant near-source grammage is encountered by CRs, or, that of a CR 
source that is farther away than in the standard leaky-box limit. 


The point to appreciate here is that, mathematically, if the simple exponential appears to be 
reasonable, data-wise that is, there is every reason to suspect that small departures and/or generaliza- 
tions in the form of Q from that of a simple exponential should still be reasonable. Or, one could even 
ask if there is a general form for Q (of which the simple exponential, and variants, are but special 
cases) that is consistent with the leaky-box ansatz. One is tempted to answer this question by going 
back to the simple one-dimensional diffusion equation and to begin relaxing the main assumptions 
that lead to the leaky-box limit as a general solution is rederived. 
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7. THE CONTACT POINT: x esc <=> D 


To recite the main conclusion of section 6 (with the question of boundary conditions aside): 
The leaky-box ansatz is the mathematical limit of the one-dimensional diffusion model for a homo- 
geneous confining volume (with the ‘right’ initial and boundary conditions) and implicit time depen- 
dence. Naturally, one needs to identify any points of contact between the two models here. The 
diffusion model, in one or more dimensions, is characterized by a diffusion coefficient (or tensor) 
D. The leaky-box limit is characterized by the escape time T esc or the escape path length x esc . Clearly 
then, a connection between the two models is manifested in a connection between D and x esc . Exam- 
ining eq. (13), one-dimensional diffusion, and equation (21), leaky-box limit, one can easily show 14 
that 




esc 


D 


(32) 


which is valid as long as h^h h 1, i.e., tube-like diffusion along the galactic magnetic held line of 
force (recall that h g is a characteristic height scale of the disc and h h for the halo). 


Equation (32) is also somewhat obvious from the second derivation of the leaky-box limit, 
i.e., the one using the potential-like function O. One can see from eqs. (29) and (30) that the constant 
must somehow be related to (a constant) D. Even though it appears absent in the literature, one 
should be able to arrive at a relation like eq. (32) using this derivation and invoking (chaotic) adia- 
batic invariance, in some simple model of (spatially) trapped charged particles. 19 


The significance of the contact point is that one can inject an energy dependence in x esc , which 
appears needed by the data, via any energy dependence in D. From a Kolmogorov-type descrip- 
tion (statistical description) of a fast charged particle’s motion in turbulent fields, 7 ' 8 ’ 20 ’ 21 the depen- 
dence of D on energy has the form 

D(E) - E~ k , (33) 

where k is a typical Kolmogorov decay exponent, characteristic of diffusive transport processes. In 
terms of rigidity R (total momentum per charge), 

D(R) - pR a , (34) 

with f5= vie and the exponent a= 1/3. With this dependence of D on R, x esc is written 

x e sc =x oP{ R 0 /R ) a ’ ( 35 ) 

with R 0 being some rigidity parameter and x 0 a constant that appears to depend on the assumed 
composition of the ISM. 22 The exponent a' seems to depend on the particular species propagated 
and whether it happens to be a primary or a secondary CR nucleus. 23 
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8. PUTTING THE BOX TOGETHER 


In this section and section 9, how one goes about implementing the leaky-box model is briefly 
sketched. From the outset, the evolution of N as a particle density distribution function has been 
discussed. What typically one measures though is the differential particle flux j (number of particles 
observed per unit area {cl A), per solid angle (cl LI), per energy interval (c/E), and per unit time (clt)). 
However, the measurement of j is equivalent to a measurement of N. What makes this true is the 
simple fact that the number of particles is a conserved quantity, i.e., 

j dE dAdQdt = N p 2 dpv dAdQdt , (36) 

where p is the nonrelativistic momentum of the particle with velocity v and mass m. Since pdp = mdE, 
eq. (36) simplifies to 

j=P 2 N. (37) 

The same can be shown for relativistic particles. 

First, reinject any energy loss/gain terms in the model for realistic numerical implementation 
of the model. Energy losses here primarily refer to those due to the ionizing effects of the charged 
CR particles as they pass through the ISM matter. Such losses are calculated using the well-known 
Bethe-Bloch expression for the stopping power W{E). Qualitatively, and at low energies, 

dF 1 

W(E) = -—ocA r . (38) 

dx \ 

At high energies, W{E) scales as lny, where yis the Lorentz factor. For ultrarelativistic particles in 
dense media, the situation is slightly more complicated, where a relativistic correction has to be 
carefully injected in the expression of W(E). 24 While W(E) is quite sensitive to the particle’s kinetic 
energy, it has no explicit dependence on the particle’s (rest) mass. This implies that particles with dif- 
ferent mass but same momentum will have similar rates of energy loss. Also, W{E) has a minimum, 
below this minimum, and because of the dependence depicted by eq. (38), the energy loss increases 
(quite rapidly) as v decreases. Above this minimum, however, the loss rate increases quite slowly. To 
a good approximation and up to two decades of E above the minimum, the loss is more or less the 
same and is of the order of a few MeV per g/cm 2 (for 1 GeV protons, for example) of matter, inde- 
pendent of the exact characteristics of the matter. 

Energy losses due to nuclear interactions — those that are responsible for species transfor- 
mation — are typically neglected; the so-called straight-ahead approximation. This is a rather good 
approximation for energies above a few GeV/nucleon. Relaxing the straight-ahead approximation 
suggests 25 that the approximation tends to introduce an error on the order of only 5%. 
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The energy loss term enters in eq. (21), the leaky-box equation, looking like 


- + ^ i [w l (E)N l (E)yq l . (39) 

Energy gain terms, e.g., diffusive acceleration, 16 stochastic acceleration, 18 or reacceleration, 26 enter 
eq. (21) in a similar manner. 


Finally, decay terms, loss as well as gain of species due to radioactive decay, also need to be 
included in eq. (21) as they tend to alter the balance of N t given q r Even though such terms are 
functions of time, and there is no explicit time dependence in eq. (21), they are taken as functions of 
energy only and a cutoff time is typically assumed, i.e., the characteristic age of CRs, on the order 
of few million years. Meaning, radioactive species that decay with half-lives longer that this time are 
considered stable here. Now, one can inject this sort of information in eq. (21) via two ways. The first 
is to write, in eq. (21), the decay term explicitly as a function of path length x as: 


clN i 

dx 


decay 


' Nj_ ^PjjNjd 

T ■ . T ■ 

1 J J J 


(40) 


where t ( - is the half-life of species i and p /; - is the branching ratio (fraction of decays of species j yield- 
ing species /'). Another way is to inject the same information about half-lives and branching ratios 
in the cross-section terms, i.e., in cr ; and appearing in eq. (2). If no explicit time dependence is 
assumed, the two ways are equivalent. 

Equation (21), in terms of the various escape, interaction, and decay path lengths, along with 
ionization energy losses, is now written as 


?;(£>- 


Hj(E) 

x' sc (£) 


M £ > , y ME 

x‘ m (E) % xf(E ) 


ME | y ME 

x?“(E) ^xf(E) 




(41) 


where the plus or minus sign of each term is obvious as the process (described by that term) either 
enhances the source term q i {+) or depletes it (-). Equation (41) is the standard leaky-box propaga- 
tion equation. 

Solving eq. (41) — for N t at a given energy, i.e., energy is treated as a parameter — proceeds in 
two equivalent ways. One can either solve for N-(x,E) in 


dN?(x,E) 

dx 


= left-handside of eq. (41) , 


(42) 
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with the source term understood as N t (x = (),E), and then integrates 


N i (E) = \*dxN°{x,E). 

Or, one solves eq. (41) for N t ( x,E ) without the escape loss term, then integrates 


( 43 ) 


r 00 

N t (E) = j o dx N- (x,E)P(x,E) , (44) 

withP(x,£) being the normalized PLD function, or Q earlier, e.g., eq. (19). The source term, 
Nj(x = 0,E) or q i (E), is typically taken to be a power law in rigidity for each species either due to 
a single shock 27 ’ 28 or multiple shocks 29 ’ 30 but with an added flattening term (power law in velocity) 
for low energies. 
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9. SAMPLE NUMERICAL IMPLEMENTATION 


In this simple — buthopefullyillustrative — example, the leaky-box model will be implemented for 
a two-species CR source made up of only protons (90%) and alpha particles ( 1 0%). At the ‘point source,’ 
these are synthesized and energized with a ‘prepropagation and spallation’ spectrum determined solely 
by the particle’s rigidity. In the ISM (made up of protons only), both the protons and the alpha particles 
suffer some energy losses but only the alphas suffer spallation. Neither suffers any decay. A simple for- 
mula will be used to account for modulation before comparing — qualtitatively only — the leaky-box- 
generated spectra to observed spectra at 1 AU. 

The two-species system of equations to be solved is then, 


and 


37V, W,(£) jV,(£) . 1 V 2 (£) 

3x q ' ( Xj esc (£) x| nt (£) x 2 ^(E) 

™ 2 = , ( E ) N 1 (E) N ^ E) i N ' (E) 

3 a * 2 ' *'«=(£) xf(£) x^ 2 (E) 


-^2(^2]. 


(45) 


(46) 


where ‘1’ refers to alpha particles and ‘2’ to protons. If the energy loss function W does not depend 
on v explicitly, which it usually does not, then the structure of this two-species system suggests the 
method of characteristics for solution. Multiply the equation for N x by W x so that it can be rewritten 
as: 


dx 


+ W X 



(47) 


where G x (x,E) is defined as W x {E)N x {x,E) and where, for the alpha particles, 


K 


1 1 

1 — r , 

esc int 

-A | A | 


(48) 


and — >°°. Using the standard method-of-characteristics solution for this linear, first-order PDE, 

gives 


A 1 (x,E) = ^^{lU 1 (Ek7 1 (E)-exp(-v-/A 1 (E)) 


W x {E')q x {E')\\-n x Xl\E') 


(49) 
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where E' = E + xW l (E). Note that the assumption is used that lim x ^ 0 N t (x,E) = nfi^E), i- 1,2, 
where n i is a constant that carries the relative abundance factor (i.e., 0.9 for protons and 0.1 for 
alphas) as well as the right dimensional conversion, without any loss of generality, the latter can be 
taken to be unity. For N 2 (x,E ), if its source term is rewritten as 

Q 2 (x, E ) = q 2 (E) + x;\ 2 /V, (x, E) , (50) 

its PDE will resemble that for N { (x,E). With G 2 defined, as before, as W 2 N 2 , 


dG 2 

dx 


+ W 2 


dG 2 

dE 


(w 2 q 2 -^‘g 2 ), 


with a solution, 


N 2 (x,E)-. 


h eg) 

W 2 (E) 


W 2 (E)Q 2 (X, E) - exp (-X / A 2 ( E)) 


x 


W^E')q 2 (E’)[\-)q\E’)+^\E’)^El 
l q 2 ( E ), 


(51) 


(52) 


where, for the protons, E' = E + xW 2 (E), and 


1 1 1 

— — 1 — 1 — r • 

3 esc int 

^2 x 2 x 2 


(53) 


Next, the various source, escape, energy loss, and interaction terms appearing in the above 
solutions for N\(x,E) and N 2 (x,E) are addressed. The source function is typically parameterized 31 
using a power law in rigidity as (see fig. 2): 


q{R) = q 0 R- 2 \ (54) 


with R = j(e + 2m 0 ^jE A/Z, m 0 being the nucleon rest mass and A and Z are the mass and charge 
number of the particles, respectively. 

The escape function is also typically parameterized 32 using a power law in rigidity as well. If 
R is expressed in units of GV, the escape path length in units of g/cm 2 can be well approximated as 
(see fig. 3): 


X cJ R > = 


26 . 7/3 

[(/ 3i?) 0 ' 5 +(/ 3i ?/ 1 . 4) -1 ' 4 


(55) 
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Figure 2. The source function, eq. (54), for protons and alphas as used 
in the sample illustration. 



Kinetic Energy (MeV/Nucleon) 

Figure 3. The escape function, eq. (55), for protons and alphas as used 
in the illustration. 


In units of MeV per g/cm 2 , the energy loss function W(E), i.e., the Bethe-Block relation, can 
be written 33 as 


W(E) = 


Afj 


V'yW r, 


max _ 


(56) 


where K= 0.307075 g/cm 2 for A = 1 g/mol, z is the charge of the traversing particle, Z the atomic 
number of the medium (1, in this case, for an ISM of protons only), A the atomic mass of the 
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medium (also 1), m e c 2 is the mass of the electron, /3 and 7 are the standard Lorentz factors for the 
traversing particle, T max is the maximum kinetic energy that can be imparted to a free electron in a 
single collision, / is the mean excitation energy, and bis a density effect term. For a traversing particle 
with mass M, 


T 

max 


2m e c 2 p 2 y 2 

1 + 2 ym e / M + m 2 / M 2 


(57) 


The density effect term, 5, can be ignored for such low density media like the ISM (see fig. 4). 



Figure 4. Energy loss by protons traversing a hydrogen gas — Bethe-Block relation, 
eq. (56). Discrete points are as simulated by the code SRIM/TRIM. 34 


For the protons-and-alphas-only leaky-box model, there are three nuclear interaction terms 
that need to be estimated. The first two are for the total reaction cross sections for CR alphas and 
protons impacting the ISM protons as a function of energy; these are needed to infer the interac- 
tion mean free paths „yj nt ( E ) and x'^iE ) . These cross sections, for kinetic energies in the range 
10 2 <E< 10 5 MeV/nucleon, are easily parameterized 2 using a simple function like: 

^total^^d = a [l + 6exp(cln/f )cos(<7 InF + c)] , (58) 

where {a,b,c,d,e} are constant parameters. For alphas, a- 100, 6=16, c = -l, d=- 2.5, and e=-2, for 
o{E) in units of mb and E in MeV/nucleon. For protons, a- 34, 6= 122, c = -l, d- 1.5, and e- 1 (see 
fig. 5). The third term is the cross section for production of protons from spallation collisions of CR 
alphas with the ISM protons, to estimate x^ 2 - This is also relatively easy to parameterize in the 
same energy range as follows (see fig. 6): 

CT pr°dC B ) = 4 1 + 6ex P (c\nE)] , (59) 
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using a- 93, 6 = -37, and c=- 1. These approximated cross sections (in mb) (typically as a function 
of energy) are then converted to mean free paths (in g/cm 2 ) using 

, r . , molar mass of the ISM [in g per mole] 10 4 

x or A [m g/cm- ] = - — — • (60) 

I Avogadro's number [in No. per mole] x o' [in cm ] ( jcj [in moj) 



Figure 5. The total reaction cross sections, eq. (58), for protons and alphas 
as used in the illustration. 



Figure 6. The spallation cross section, eq. (59), for alphas on ISM protons 
to produce secondary CR protons as used in this illustration. 
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Next, the solutions need to be averaged, which are now functions of both distance or path 
length (in g/cm 2 ) and energy using an assumed PLD. In the leaky-box approximation, this distribu- 
tion, eq. (19), when normalized becomes 


Q(x) = x es c ex p(-A‘ / x esc ) , 


(61) 


where v esc is the mean escape path length (typically <10 g/cm 2 ). Hence, the arriving solutions (but 
still unmodulated yet by the solar wind) become (z = 1,2), 

(£) = J 0 °° N i ( X ’ E Ksc ex P (~ x 1 *esc ) dx • ( 62 ) 

Solar modulation, for our purposes, can be illustrated using a potential 35 and an energy loss as, 


E +2m n E - 

NAE) = ~ 5 NAE + AE), 

1 (E + AE) 2 +2m Q (E + AE) 1 


(63) 


where AE = Ze<f>/A, O being the potential (in V; typically about 400 MV for solar minimum condi- 
tions and about 800 or higher for solar maximum conditions) and e the electron charge. 

The ‘arriving’ number densities at 1 AU for protons and alpha, N 2 (E) and N^E), are 
shown in figure 7, where the averaging over the path length, eq. (62), has been done numerically. 
Here, no atmospheric effects are taken into account. In figure 8, the spallation of alphas’ contribu- 
tion to the arriving protons are highlighted. This easily is calculated using 100 x^N 2 (E)l N 2 (E)), 
where N' 2 (E) is a solution to the proton equation, eq. (51), but with the spallation term ‘turned off’ 
by letting Figure 9 shows the ratio of the arriving alphas to the arriving protons, i.e., 

N\(E) and N 2 (E ) , which includes the secondary proton component. 



Figure 7. The arriving number densities for protons and alphas, eq. (63), as predicted 
by this sample illustration using x esc - 1 g/cm 2 and 0=400 MV. 
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Figure 8. Ratio of the secondary protons (produced by the spallation of the primary CR 
alphas in the ISM) to the primary CR protons. 
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Figure 9. Ratio of arriving alphas to protons (including secondary protons). 


To put these sample calculations in perspective, figure 10 shows actual observed BESS 36 data 
for the flux of CR protons and alphas; their ratios are shown in figure 11. BESS balloon-borne spec- 
trometer was flown from Manitoba, Canada, in July 1998 (i.e., under solar minimum conditions) at 
a height of ~37 km, where residual atmospheric effects are equivalent to <5 g/cm 2 . For this comparison, 
the measured BESS flux has been multiplied by a constant (= 2x 10 -8 ) so as to scale it to the calculated 
number density. The results from this simple illustration are not meant to be a model prediction designed 
to be compared directly to data. Only qualitative comparisons are intended and useful here. Already, 
a number of general observations can be made: 
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Figure 11. Ratio of alphas to protons as observed by BESS. 36 
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(1) The general shape of the protons’ and alphas’ arriving spectra appear straightforward to 
reproduce assuming simple rigidity dependent escape and source functions. 

(2) Alphas do not appear to be a major contributor to secondary protons. 

(3) Spallation of CR nuclei heavier than alphas must contribute to both protons and alphas. 

(4) Low energy (less than a few GeV’s per nucleon) behavior of GCRs is more difficult to 
model than at higher energies, where solar modulation effects and nuclear (total reaction and spall- 
ation) cross sections appear to subside. 
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(5) Propagation effects (nuclear chemistry of CRs) appear to be well captured by simple mod 
els like the leaky-box model mainly because the energy dependence of the relevant nuclear param 
eters is rather weak as well as being well understood and modeled. 
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10. REMARKS 


In ending this brief review focusing on the practical side of modeling the galactic propagation 
of cosmic rays, the following are three general closing remarks: 

(1) The intended reader is a newcomer to the held as well as that practitioner who simply 
wants to delve a bit deeper into this immensely interesting area. At various points in this review, I 
suspect both wondered why I had focused and picked areas to elaborate on and paid only scant atten- 
tion to others. I have no strong general justification to this (accurate) observation except to say the 
following. I had ventured my own study into cosmic rays from a nuclear physics background (and 
predisposition). Learning the eclectic new subject, I had to build my knowledge base and formulate 
my own understanding from what seemed to me like a mosaic of specialities, styles, backgrounds, 
rigors, and aims. It was ‘natural’ for a nuclear physicist to gravitate towards the nuclear aspects of 
cosmic rays. In retrospect, I believe, it was equally ‘natural’ to also look for something ‘different’ to 
do. Galactic propagation aspects of cosmic rays are related (cf. nuclear reactions) yet are somewhat 
different (cf. kinetic and transport formulations). 

(2) I was attracted by the Ginzburg- Syrovatskii elegant formulation of the problem from 
early on. I had found their derivation and exposition of the kinetic theory of this problem (‘the 
quantitive theory of the origin of cosmic rays,’ as they had auspiciously termed it) an area I wanted 
to get into. Their theory survives to this day. I strongly recommend that the interested reader pick 
up and study their classic book. The limit of the Ginzburg- Syrovatskii propagation theory, the 
leaky-box model, is now a standard tool for both cosmic-ray data analysis and modeling as well as 
a more or less complete model of the nuclear chemistry of cosmic rays. In other words, the leaky-box 
model, albeit a limit to or a crude approximation of ‘the quantitative theory of the origin of cosmic 
rays’ is essentially a complete framework for the nuclear chemistry of cosmic rays. 

(3) Detailed modeling — with an eye on data comparison and such — of the nuclear chemistry 
of cosmic rays using the leaky-box framework is actually quite straightforward as attempted to be 
demonstrated. ‘Professional’ versions of this leaky-box demonstration (e.g., see ref. 37) are funda- 
mentally similar to this simple one. 

The first 100 years since the discovery of cosmic rays, and especially since the coming of 
the space age, has seen unimaginable progress, refinement, and evolution of understanding of this 
intriguing phenomenon. It is quite likely, however, that the nuclear chemistry aspects of cosmic rays 
may have matured to a point where similar and equally fantastic progress cannot be foretold — in yet 
another sigmoidal display of change in nature. 
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